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Abstract 


In order to decrease overall computational time requirements of a spatially-marching 
Parabolized Navier-Stokes finite-difference computer code when applied to turbulent fluid 
flow, a wall-function methodology, originally proposed by R. Barnwell, was implemented. 
This numerical effort increases computational speed and calculates reasonably accurate 
wall shear stress spatial distributions and boundary-layer profiles. Since the wall shear 
stress is analytically determined from the wall-function model, the computational grid near 
the wall is not required to spatially resolve the laminar-viscous sub-layer. Consequently, 
a substantially increased computational integration step size is achieved resulting in a 
considerable decrease in net computational time. This wall-function technique is demon- 
strated for adiabatic flat plate test cases from Mach-2 to Mach-8. These test cases are 
analytically verified employing: (1) Eckert reference method solutions, (2) experimen- 
tal turbulent boundary-layer data of Mabey, and (3) finite-difference computational code 
solutions with fully resolved laminar-viscous sub-layers. Additionally, results have been 
obtained for two pressure-gradient cases: (1) an adiabatic expansion corner and (2) an 
adiabatic compression corner. 
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1. INTRODUCTION 


1.1 Research Topic 

To accurately resolve both the viscous drag and heat transfer effects of a compress- 
ible turbulent flowfield, standard computational fluid dynamics (CFD) techniques require 
prohibitive amounts of computational time for problems of engineering interest. Numer- 
ous researchers have demonstrated that to accurately calculate a turbulent boundary layer 
employing finite-difference computational techniques at least one grid point must reside 
in the laminar-viscous sub-layer (i.e. the inner portion of a turbulent boundary layer) [1], 
Hence, for turbulent boundary-layer calculations, the first grid point off the wall must be 
at a y + <l-2 (denoted as the fully-gridded CFD case throughout the text). The parameter 
y + is the transformed coordinate of the normal-wall coordinate, y, defined by, 



where r w is the wall shear stress, p w is the wall density, and /r w is the wall viscosity. 
When this y + constraint is applied to a uniform grid-generating scheme, typically a 
computational grid is generated with hundreds of grid points in the boundary layer, 
which ultimately requires a very small integration step size between solution planes. 
Numerous investigations have addressed the grid spacing problem with grid stretching 
algorithms being one of the most commonly employed solution methodologies [1], Grid- 
stretching algorithms yield non-uniform grids with grid points clustered in high gradient 
regions. Unfortunately, this approach only moderately influences the computational 
time expenditures necessary to derive meaningful engineering calculations. Hence, wall- 
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function methods are utilized to substantially reduce the computational time necessary to 
generate solutions for turbulent flowfields. 

1.2 Wall-Function Methods 

In general, wall-function methods calculate an analytic inner-flowfield solution and 
patch this inner solution to a numerically-generated outer-flowfield solution at a location 
denoted as the match point. At the match point, the shear stress, velocity, and turbulent 
viscosity (or equivalent quantities) are matched, implying that the computations of these 
quantities using a wall-function method are used as boundary conditions for CFD codes. 
Typically, the match point is in the logarithmic region of a turbulent boundary layer, 
between a y + of 40 and a y + of 400, seen for example in Figure 1 for incompressible 
flow, taken from reference [1]. 



y + 


Figure 1: Law-of-the-Wall Velocity Profile. 
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The non-dimensional velocity, u + , is defined as. 


where u* is the shear-stress velocity (commonly denoted as u r ) and Re$ is the Reynolds 
number based on momentum thickness. Researchers have developed methods to ana- 
lytically calculate the inner region of the boundary layer without having to explicidy 
solve the equations of motion in that region. In short, the law-of-the-wall is the basis 
for wall-function methods. Numerous methods require that the user specify the location 
of the match point to lie within the inner portion of the boundary layer at each stream- 
wise location. These approaches do not easily allow for the optimum placement of the 
match point [2], 


1.3 Defect Wall-Function Method 

The Barnwell and Wahls [3, 4, 5, 2] wall-function method was developed for analysis 
of both incompressible and compressible flows, adiabatic and non-adiabatic flows, and 
zero pressure-gradient and non-zero pressure-gradient cases. This method differs from 
previous wall-function methods in many aspects. The Barnwell and Wahls method uses 
analytic-velocity functions, the law-of-the-wall, the law-of-the-wake, and the defect- 
stream function, to calculate wall shear stress and slip-wall boundary conditions consistent 
with the inner and outer analytic-flowfield solutions. Since the wall shear stress is 
analytically calculated, the first grid point off the wall does not lie within the laminar- 
viscous sub-layer, thereby eliminating the prohibitive grid spacing constraint. The slip- 
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wall boundary conditions are employed by numerical codes to generate solutions by 
computing the flowfield all the way to the wall and thus, there is no need to patch 
an analytic-inner solution to an outer solution generated numerically. Consequently, 
the computed inner layer is non-physical and an extension of the outer layer to the 
wall. This eliminates the need for an inner turbulent eddy-viscosity model [2], hence, 
the outer turbulent eddy-viscosity model generates viscosity values at all numerical grid 
points. Previous methods employed the inner-layer analytic-velocity functions as the 
inner-flowfield solution and patched a numerically calculated outer-layer solution to the 
analytically calculated inner-layer solution at the match point consistent with the boundary 
conditions at the match point. 

The grid point where the inner and outer layers meet is denoted as the match point. 
One of the advantages of the Barnwell and Wahls wall-function method is specifically 
related to the calculation of the match point location. Specifically, this wall-function 
method calculates the location of the match point, whereas numerous other methods 
required the user to specify the match point location (note that the match point is a function 
of streamwise location). The Barnwell and Wahls method self adjusts the location of the 
match point at each streamwise position and therefore is easier to implement, since no 
user specified information about the match point is required a priori. Additionally, the 
Barnwell and Wahls [2] match point is forced to be at the optimum location (i.e. the outer 
edge of the inner layer) and thus potentially allows for larger grid spacing as compared 
to those of previous wall-function methods. 

The analytic functions employed by this wall-function method are: (1) the law-of- 
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the-wall and law-of-the-wake for the inner portion of the boundary layer, and (2) the 
defect-stream function for the outer portion of the boundary layer. These functions 
are discussed in detail in Section 2.2. Another advantage of the Barnwell and Wahls 
method is the use of the law-of-the-wake as part of the inner-layer velocity function. 
The law-of-the-wake extends the effective region of the law-of-the-wall and also allows 
the streamwise pressure gradient to influence the inner region of the boundary layer [2], 
The law-of-the-wake has been used in other research to describe the velocity profile for 
the outer region of the boundary layer [6], but in this development it is only used for 
the inner region. In short, the Barnwell and Wahls wall-function method uses only the 
defect-stream function to evaluate the outer-layer analytic-velocity profile. 

The defect-stream function is based upon studies by Clauser [7] of equilibrium 
turbulent boundary layers. Clauser defined a boundary layer to be in equilibrium if 
the following condition is satisfied, 

< 5 ; d P 

— = constant , (3) 

t w dx 

where 8* is a boundary- layer displacement thickness parameter and is the streamwise 
pressure gradient. This condition, if satisfied, represents a balance between the pressure 
forces and the shear forces in a turbulent boundary layer [8] and is assumed to be valid 
for all turbulent cases analyzed in this research. The zero pressure-gradient case is a 
special case of equilibrium boundary-layer flow. 

In the Barnwell and Wahls wall-function method a slip (non-zero) streamwise veloc- 
ity, consistent with the analytically calculated wall shear stress, is imposed at the wall to 
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permit integration to the wall. In contrast, standard CFD codes explicitly set the stream- 

wise velocity on solid surfaces to zero. To determine the slip-wall velocity, the velocity 

gradient at the wall is calculated from the definition of shear stress, 

du 
dy 

where the laminar viscosity, /*/, is calculated employing Sutherland’s law [9] and the tur- 
bulent eddy viscosity, /q, is calculated using the Baldwin-Lomax turbulence model [10]. 
A first order finite-difference approximation of the velocity gradient is used to define the 
consistent slip streamwise velocity at the wall, u s , 

l ")\l 

u a = u(l) = u(2) — [y(2) — j/(l)] , (5) 

where u(l) is the streamwise velocity at the wall, u(2) is the streamwise velocity at the 
first grid point off the wall, y(l) is the y-location of the wall, and y(2) is the y-location 
of the first grid point off the wall. The slip-wall velocity is subsequendy utilized as a 
boundary condition for a CFD code to numerically calculate the entire flowfield (i.e. a 
CFD code integrates the entire distance to the wall). Furthermore, a slip-wall density and 
slip-wall temperature are calculated consistent with the slip-wall velocity. 

This research effort investigated the application of the Barnwell and Wahls wall- 
function methodology to reduce the time requirements of a Parabolized Navier-Stokes 
(PNS) CFD code, developed by Korte [9]. The code uses an explicit, upwind, space- 
marching finite-difference scheme to eliminate time as a variable and permits the use of 
a non-iterative or single-pass technique to resolve the flowfield. The PNS equations are 
commonly utilized (when relevant), instead of the full set of Navier-Stokes equations, 
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since less computer memory and less computational time are needed to generate solutions. 
The PNS equations have been shown by researchers to accurately calculate flowfields 
within reasonable time constraints [ 1 ]. 

1.4 Computational Fluid Dynamics Code 

The finite-difference code used in this research effort solves the Parabolized Navier- 
Stokes (PNS) equations with an explicit, upwind space-marching scheme [9]. The PNS 
equations are derived from the full set of unsteady Navier-Stokes equations by neglecting 
the unsteady terms, neglecting the stress and heat flux terms with respect to the streamwise 
direction, and neglecting a fraction of the subsonic streamwise pressure gradient by 
employing Vigneron’s coefficient, u>. Vigneron’s coefficient has a valid range of 0.0-1.0 
and is applied with a safety factor, a, in the following form, 

U> = min (1, <7u>) , (6) 

where 

(Mf > 1) 

2 , (7) 

tfflj 

where a typical value of a has been taken to be around 0.75 for this research, is the 
axial Mach number, and 7 is the ratio of specific heats. The PNS equations are a mixed 
set of hyperbolic-parabolic differential equations assuming that the inviscid portion of the 
flow is supersonic and the streamwise velocity is positive. The latter constraint demands 
that the flow be attached at all streamwise locations (i.e. streamwise flow separation is 
not permitted) [ 1 ]. 
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The PNS equations, in general, are applicable to two or three-dimensional, steady, 
supersonic, viscous flowfields without streamwise separation. The advantage of the 
implementation of the PNS equations compared to the full set of Navier-Stokes equations 
is that the solution is obtained with an efficient space-marching method, producing faster 
execution times and using less computer memory. 

The non-dimensional form of the two-dimensional PNS equations, used in the finite- 
difference code developed by Korte [ 9 ], are presented in the transformed coordinate 
system, £-77. The transformed coordinate system was developed to handle complex 


geometries. The governing PNS equation is: 


where 


+(¥ +2 £ 




(Vv\ 

\ J ' T) 


,*'(n 


£=£(**) 77 = T](x*,y*) , 


E = Ei- E v F = F, - F v , 


E, = E* + P , 


p* 

p* U* V* 
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(I-w)jj* 

0 
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( 15 ) 


F t = 


p*V* ^ 

p*v*u* 
p m v*v*+p m 
(e;+p*)u* ^ 


*> 


F v = 


u*r; y +v*T* y -ql ) 


( 16 ) 


* 

e, = 


e* + 


u 


+ u* 


( 17 ) 


where e is the internal energy, e t is the total energy, J is the Jacobian of the transformation, 
p is the pressure, q is the heat flux, u and v are the velocity components, p is the density, r 
is the shear stress, and the superscript * implies a non-dimensional quantity. The subscript 
“i” denotes an inviscid parameter and the subscript “v” denotes a viscous parameter. 


The variables in equations (9) through (17) have been non-dimensionalized with the 
following relations, 


.* _ x 
L 


y = 


p* = 


U oo 

p 

Poo Use 


V = 


U o 


Poo 
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( 18 ) 


T* = 


p OoU( 


o 

oo 


e* = 


U 2 

w oo 


where oo denotes freestream conditions, L is the characteristic length, and T is the 
temperature. 


The parabolized forms of the shear stress, r, and heat flux, q, terms in the transformed 
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coordinates for the above equations are: 


= 


T jcy 


'[2(v*«?|) - (Vy v *,)} 


* 

<h 


yy 


2/x* 

3 Re L 

[('/*<) - M)] 

-/** 
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( 19 ) 


V 

3 rile 


% = 


h - l)Ai:i,Rc L Pr 




where 


/<* = N + /T 


( 20 ) 


and ReL is the Reynolds number based on the characteristic length and Pr is the Prandtl 
number. The non-dimensional laminar viscosity is calculated using Sutherland’s equation 
shown below, 


f L i 


3 

T *~2 


/IT T r(: f \ 
\T* T T~f) 


(21) 


where 

110.4 K 

T re f = 7p 

r oo 


( 22 ) 


The calculation of the turbulent viscosity is discussed in Section 2.1. 

One of the important features of this finite-difference code is the non-uniform grid 
capability. The grid points are clustered near the wall to ensure adequate resolution of 
the laminar-viscous sub-layer where the gradients between grid points are large, and a 
sparser grid is used farther from the wall where the gradients are smaller. The grid point 
locations are generated employing a typical Robert’s stretching function [1], which has 
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been modified by Korte to have the following form. 


y(j) ~ Drain H" ( Umax ~ 2/min)^l — sfact + 2 — — ^ j , (23) 

where ymin ^nd y m ax <trc the y-locations of the boundary of the computational domain. 
One of the parameters to control the stretching, sden, is a function of the grid point 
number, j, and the stretching factor, sfact, given as, 

, , f s fact + 1\ 

sien = i + {jjzr^T ) • <24 > 


where nj is the total number of grid points in the field. The stretching factor controls 
the clustering of the numerical grid near the wall. Typical values for sfact range from 
1.01 to 1.000001 depending on the grid resolution required. Examples of the effect of 
the stretching factor on grid point location using the same total number of grid points 
are presented in Figure 2. 



y/L 

Figure 2: Grid Point Distribution for Varying Stretching Factor. 

If the stretching factor is changed from 1.01 to 1.0001, it effectively puts 21 grid points 


it 



close to the wall as compared to only 2 for the 1.01 case for the same flowfield height. 
The stretching factor is calculated at each streamwise location to maintain boundary- 
layer resolution in the laminar-viscous sub-layer as the code marches downstream for the 
fully-gridded CFD computations. The y + location of the first grid point off the wall (user 
specified prior to compilation of the code) typically has a value of y + <l— 2 to resolve 
the laminar-viscous sub-layer for turbulent flows. 

The second-order accurate, two-stage (i.e. predictor-corrector), explicit, upwind 
scheme used by Korte to solve the PNS equations is: 


Stage 1: 


= -(V;\ - fi _j+ 


(n + > - n) + VCD” + (fy) (n? - o;- 1 ) 




T p -T p ) + 


(n -n-) + < cc£ >5 + 


~p - !I± e> + 1 JjLf l 

1 Vj — J J-'Vj ^ J v ) ’ 


Q — w 

H — 0 


and the superscript n represents the values at a known flowfield plane, p represents the 
predictor stage values, and n+1 represents the unknown flowfield plane to be determined. 
Basically, Stage 1 calculates the p values using the n values, then the p and n values 


are used to calculate the n+1 values in Stage 2. The prime denotes viscous stress and 


12 



heat fluxes to be differentiated with respect to the ^-direction. The GCL parameter is 
the Geometric Conservation Law term defined as. 


{GCL)] = Ef 



0 GCL )J = e; p 




(31) 


where 

E' = Ei - E' v 
F' = Fi - F' v . 

The upwind flux approximations are obtained using Roe’s flux-difference splitting method 
by either splitting the flux vectors or flux differences based on the sign of the eigenvalue 
(or wave speed). Roe’s method has been modified for the PNS equations by Korte [9], 
Note that there was no modification to the integration scheme developed by Korte for 
this research. 


1.5 Implementation of the Defect Wall-Function Method 

Explicit finite-difference solution methods have a more stringent stability constraint 
(i.e. the CFL number) than implicit methods, restricting the step size; therefore a wall- 
function method is implemented to increase the integration step size. The implementation 
of the wall-function method dramatically reduces the computational time requirements 
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of the CFD code, since the grid spacing is increased by eliminating the need for the y + 
constraint, and consequently produces a much larger integration step size. 

The Barnwell and Wahls wall-function method is implemented in the Korte finite- 
difference PNS CFD code to analytically calculate the wall shear stress and to determine 
consistent slip-wall boundary conditions for utilization in the PNS CFD code. The 
finite-difference code with no-slip boundary conditions is employed for the first few 
computational planes (user specified) in order to generate initial data for the application 
of the wall-function method. As the PNS CFD codes marches, the previous computational 
plane is used as the initial solution for the predictor stage (and the solution for the predictor 
stage is the initial solution for the corrector stage), consistent with the computational 
algorithm. The wall-function methodology in the modified finite-difference code is used 
exclusively to advance the spatial marching procedure. 

The Barnwell and Wahls wall-function method was originally developed using the 
Clauser turbulence model (characterizing the outer portion of the boundary layer). The 
PNS CFD code uses a Baldwin-Lomax turbulence model for the outer portion of the 
boundary layer. Hence, the Barnwell and Wahls wall-function method was modified 
in this research to utilize the Baldwin-Lomax turbulence model. This is discussed in 
Section 2.1. In this research, the freestream conditions (imposed upon the grid point at 
the outer edge of the computational domain and denoted with the subscript oo) are used 
to approximate the boundary-layer edge conditions. Note that this assumption is exact, 
except for the non-zero pressure gradient cases which are detailed in Section 4.3. 

The first grid point off the wall is placed at some percentage of the match point 
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location, as specified by the user and the remaining grid point locations are generated 
according to the original stretching algorithm developed by Korte [9]. A discussion of 
this is in Section 3.2. Since the computed flow between the match point and the wall is 
just an extension of the outer layer to the wall in the Barnwell and Wahls wall-function 
method, the outer-layer turbulence model is applied all the way to the wall. 

An “analytical grid”, a collection of discrete y-locations (not related to the numerical 
grid), is generated within the boundary layer and allows the calculation of the vorticity 
distribution using the analytic-velocity formulations (the vorticity is required by the 
Baldwin-Lomax turbulence model). Vorticity is a function of the y-location and is 
approximated by, 

du 

w « -p , (32) 

oy 

(in this research) for two-dimensional flat plate flows. The analytic-grid has adequate 
resolution in the boundary layer to calculate a reasonably accurate vorticity distribution, 
thus allowing for the utilization of a sparser numerical grid. This is discussed in 
Section 3.3. 

The integration scheme of PNS CFD code was verified using laminar flat plate 
test cases. The results of these test cases compared well to the theoretical analysis 
of Crocco [11] and confirm that the code is functioning properly for laminar flow. These 
results are presented in Section 4.1. The CFD code was modified to incorporate the 
wall-function methodology and then tested at several different Mach numbers employing 
turbulent flow conditions using both zero and non-zero streamwise pressure gradient 
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cases. The modified code produces results more efficiently than the fully-gridded CFD 
code. An order of magnitude increase in speed was obtained for an adiabatic Mach-2 flat 
plate case with only a 15% difference in the calculation of the wall shear stress (compared 
to the fully-gridded CFD solution). The modified code has also been shown to calculate 
with moderate accuracy, the boundary-layer profiles (velocity, temperature, and density). 
Additionally, the conservation of mass, momentum, and energy was checked and found to 
be reasonable for both the fully-gridded CFD case and the wall-function method. These 
results are presented in Section 4.2. 

The slip-wall velocity and density, the latter based on the empirical formula of 
Crocco [3], numerically produces a small non-physical streamwise mass flux at the 
wall. The implications of this inherent property of the Barnwell and Wahls wall-function 
formulation were not fully addressed in this research, but rather the implementation of 
this methodology in a practical computational scheme. 

The implementation of the Barnwell and Wahls wall-function method has been 
proposed to relax the grid resolution constraint for analysis of turbulent fluid flows within 
the Korte PNS CFD code. The concept of wall-functions has been introduced as well as 
the basic ideas of the Barnwell and Wahls wall-function method. The PNS CFD space- 
marching code, developed by Korte, has also been introduced in order to understand some 
of the concepts dealt with in this research, such as the predictor-corrector integration 
stages and the computational gridding scheme. Also included was a discussion of the 
modifications required to apply the wall-function theory to the PNS CFD code. The 
subsequent text presents a detailed description of the Barnwell and Wahls wall-function 
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method and the modifications to the method required for adaptation to the Korte PNS 
CFD code. Also included is a discussion of the modifications to the CFD code required 
in order to implement the wall-function method. Results from the validation of the code 
for laminar and turbulent flows are presented as well as the results from the application 
of the wall-function method and the comparisons with the fully-gridded CFD code. 
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2. DEFECT WALL-FUNCTION METHOD 

The Barnwell and Wahls wall-function method has distinct advantages compared 
with other wall-function methods. One pertinent advantage is the development of an 
equation to specify the location of the match point at each streamwise location. Another 
advantage is the use of analytic-velocity functions enabling the direct calculation of the 
wall shear stress and the corresponding slip-boundary conditions employed by the PNS 
CFD code to numerically calculate the entire flowfield to the wall. This method ultimately 
relaxes the grid resolution constraint (compared with fully-gridded numerical schemes), 
thus allowing larger integration step sizes to be employed. A discussion of the relevant 
Barnwell and Wahls wall-function theory is presented in this section. 

2.1 Tiirbulence Models and Match Point Equation 

Barnwell and Wahls [3] derived the match point equation by equating the inner 
and outer-layer turbulent eddy-viscosity models. The derivation of this equation used 
a Prandtl/Van Driest turbulent eddy-viscosity model for the inner region [10, 12] and a 
Clauser turbulent eddy-viscosity model for the outer region [3]. However, the original 
PNS CFD code employs a Baldwin-Lomax [10] turbulent eddy-viscosity model for the 
outer layer, thus the match point equation was rederived for this research using this 
turbulence model. Presented in the following text is a discussion of all the relevant 
turbulence models. 
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a) PrandtI/Van Driest Inner-Layer Turbulence Model 


Prandtl derived a mixing-length formulation based on a simple physical model of 
the turbulent shear stress. 


Tt = 


—p u'v' 



du du 
dy dy 


(33) 


where u' and v' are the time-averaged velocity fluctuations and / is the mixing-length 
parameter [12]. The mixing length is analogous to the mean free path between molecules 
of a gas; Schetz has said, “[the mixing length] is taken as some effective interaction 
distance, except that it is between eddies rather than molecules” [12]. Van Driest derived 
a mixing-length model for the inner portion of the boundary layer (the laminar- viscous 
sub-layer, the buffer zone, and the law-of-the-wall region) and is defined as, 


/ = 


Ky 


1 — exp 



(34) 


where the constant A + is equal to 26 as suggested by reference [10] and the von Karman 
constant, k, has a value of 0.41, as suggested by reference [4]. A predicted velocity 
profile for a turbulent boundary layer using the Van Driest mixing-length turbulence 
model for the inner region is presented in Figure 3, taken from reference [12]. 
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Figure 3: Inner-Layer Velocity Computed Using Van Driest Mixing-Length Model. 


The Prandtl/Van Driest mixing-length turbulence model [10] defines the turbulent inner- 
layer eddy viscosity as, 


Vt = P? M 


(35) 


where the magnitude of the vorticity is: 


M 



dv\ 

dx) 


+ 


/ dv dw \ 2 / dw du \ 2 
\dz dy ) + \ dx dz ) 


(36) 


for three-dimensional flows. 


b) Clauser Outer-Layer Turbulence Model 

Clauser developed an outer-layer turbulence model based on an eddy-viscosity model 
derived from a generalized defect-law formulation [12] presented in Figure 4, detailing 
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several transformed velocity profiles (presented in Figure 5). 


ti Vu'(0)/(Uoo(n)5*) 


0 12 3 4 



The turbulent transport coefficient for the model is assumed to be constant across the 
outer region, 


/<« = M*) + f(y) , (37) 

and the turbulent velocity profile typically intersects the wall at a non-zero value [12] 
as shown in Figure 5. 
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Figure 5: Clauser Velocity Profiles with Non-Zero Wall Velocities. 

Clauser proposed employing a pseudo-laminar (i.e. constant turbulent eddy viscosity) 
outer boundary layer model and from this derived an equation of the same form as 
Blasius’ laminar flat plate solution, which (when properly transformed) collapses all the 
data sets onto a single curve, as previously presented in Figure 4. 

Also, on dimensional grounds, Clauser proposed a turbulence model, 

/i< ex density * velocity * length , (38) 

where the characteristic velocity was chosen as the shear-stress velocity and the char- 
acteristic length was chosen as the integral thickness, A ,• [12], based on the defect-law 
formulation, 

1 Uoo - u(y) d (y\ 

u* \6 ) 




A, 

6 
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(39) 



where U 0 0 is the freestream velocity and 8 is the boundary-layer thickness. Hence, the 
turbulence model is: 


fit = kpu* A, , 


(40) 


where 



(41) 


(42) 


The outer viscosity model developed by Clauser is presented in equation (43), where the 
Clauser constant, k, is typically set to 0.0168 as suggested by reference [10], 


pi — k p Ij og 8 t 


(43) 


where 8* is the incompressible boundary-layer displacement thickness parameter, defined 


as. 



(44) 


This turbulence model was utilized in the original development of the Barnwell and 
Wahls wall-function method. 


c) Baldwin-Lomax Outer-Layer Turbulence Model 

The finite-difference computational code developed by Korte employs the Baldwin- 
Lomax [10] turbulence model for the outer-layer turbulent eddy viscosity. The Baldwin- 
Lomax model, 

fit k p Ocj> F wake F kit b ( !J ) ) (45) 
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replaces the Clauser model in the outer layer for this research. Both models have a 
similar form and yield identical //( values if is replaced by C cpF wa ic e Fkieb{y ) • The 

parameters used in the Baldwin-Lomax turbulence model are: 



where the variable ymax is the y value where F(y) is a maximum (i.e. F max )> 


u dif 


(\/ U- + V- + u.’-'j - (x/V 2 + V 2 + u» 2 ^ 


(48) 


and 


Fuib(y) 


1 + 5.5 


( b y \ b 

'Umax J 


(49) 


The Fkieb parameter is the Klebanoff intermittency factor characterized as the fraction 
of time that a flow is locally turbulent. Thus, this factor causes the Baldwin-Lomax 
model to yield turbulent viscosity values tending toward zero far from the wall. The 
coefficients employed in the Baldwin-Lomax turbulence model are listed in Table 1 as 
recommended by reference [10]. 


Table 1: Coefficients for the Baldwin-Lomax Turbulence Model. 


A + 

Cep 

Cklcb 

C\vk 

k 

26.0 

1.6 

0.3 

0.25 

0.0168 
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d) Derivation of Match Point Equation 


The match point equation originally developed by Barnwell [3] is derived by equating 
the Prandtl/Van Driest inner-layer turbulence model, equation (35), with the Clauser 
outer-layer turbulence model, equation (43), 

kpUooS* = pi 2 M . (50) 


The match point equation is used to determine the transition y-location between the inner 
and outer-layer models. The derivation of the Barnwell match point equation follows 
in the subsequent text. For two-dimensional flat plates, the magnitude of the vorticity 
is approximated by. 


jo; | ~ 


dii 

Oy 


(51) 


and is calculated using the analytic-velocity profile from the law-of-the-wall (g(y + ), 
equation (76)) and law-of-the-wake (h(y + ), equation (77)) formulation. 


u = u*[tf(2/ + ) + /«(y + )] 


(52) 


The mixing-length parameter, /, is approximated by, «y, where k is the von Karman 
constant, hence, 


kpU^S* = p{kij ) 2 ^-{u*[(/(y + ) + h(y + )]} . (53) 

The shear-stress velocity is a constant with respect to y, thus the equation reduces to 
the following form. 
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( 54 ) 


kpU ooS: 



Ojj{y+) dh{y+) 
dy dy 


The law-of-the-wall and law-of-the-wake advocated by Barnwell and Wahls are substi- 
tuted into the equation, yielding, 



(note: the wall-function parameters, R, W, t>*, and // are discussed later, in detail, 
following the match point equation derivation). The transformation from y + coordinates 
to y coordinates is governed by, 


y_ _ Pw 

A P W^OO^V 


( 56 ) 


where A is the boundary-layer thickness parameter. After taking the spatial derivatives, 
the following equation results, 


kpUooSt 



— In j/ + + 6 

K 


V 

h 

Ky 


/?( sin 


v{ — In y + + 6 

K- 


V 

KIJ 


+ 


12 y 
— W— 

vv a "> 
K A- 


) 


( 57 ) 
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The transformation from y coordinates to ?/ coordinates is governed by, 



( 58 ) 


thus, 


kUppSt 

HU* A 


f- 

' 

l" 

cos 


v ( — In y + + b 

K 


+ 


R\ sin 


v ( — In y + -f b 


Barnwell defines a parameter, u> w r, such that. 


+ [12 Wif] 


(59) 




(60) 


where 6* is the density-weighted velocity thickness parameter (as defined by 
Barnwell [3]), 


= 


JO Poo 

An alternate form for 6* is presented as, 

K = — 

II 

Uoo 

and is substituted into equation (60) to yield, 




U Q 


(61) 


(62) 


LV 


wf 


Uoo6*i 
u* A 


(63) 


Further substitution employing equation (59) yields the Barnwell and Wahls match point 
equation: 
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( 64 ) 


121V)/ 3 + y 



+ 


R \ sin 


v ( — In y + -f b 

K 



W wf k 

K 


This equation has three roots, but only one root is both positive and real. Solving this 
equation yields the match point location in the ?/ coordinate system. In y + space, the 
match point location is designated as y m + , with typical values being on the order of 
10-1000 for the test cases of this investigation. 

The wall-function parameters used in the match point equation are: 



( 65 ) 

( 66 ) 

( 67 ) 

(68) 
( 69 ) 


where /? is the compressible pressure-gradient parameter defined as. 



1 + — - + (1 — r) 
Poo 



The incompressible pressure-gradient parameter, /?, is: 


= 


K 

p w u * 2 dx 


( 70 ) 


( 71 ) 


28 



The variable T aw is the adiabatic wall temperature, 


Tnw = T a 




( 72 ) 


(assuming a constant ratio of specific heats, 7) and r is the recovery factor given by the 
empirical turbulent correlation, 


r = Pr * 


(73) 


The current research utilized the Baldwin-Lomax turbulence model instead of the 
Clauser turbulence model, hence the match point equation was rederived. The derivation 
of the match point equation using the Baldwin-Lomax turbulence model is similar to 
the derivation of the match point equation using the Clauser turbulence model. The 
terms from the Clauser turbulence model, U^S*, must be replaced by the terms from the 
Baldwin-Lomax turbulence model, C Cj> F wakc F kU:h {y), in equation (59). The Fkieb term 
is near 1.0 around the match point and thus is neglected. The governing match point 

equation utilizing the Baldwin-Lomax turbulence model is: 

kC cp^wake 


Ml 


‘A 



cos 

i/( — lnj/ + + 6) 

l 

. 

\« )\ 


+ 


( 74 ) 


1 


v I — In y + -f b 

K 


R (sin 

and upon substitution and rearrangement is: 


+ 


[1211V]} , 



This is the match point equation employed by the modified PNS CFD code to determine 
the match point location. 

Some examples of the match point distribution for the flat plate test cases are shown 
in the following figure for a range of Mach and Reynolds numbers. 



Figure 6: Match Point Distribution for the Flat Plate Test Cases. 


2.2 Analytic-Velocity Functions 

The analytic functions employed for the inner layer are the law-of-the-wall, equation 
(76), and the law-of-the-wake, equation (77). The analytic outer layer is governed by 
the defect-stream function, equation (80), a measure of the deficiency of the outer-layer 
velocity (compared to the freestream velocity). The law-of-the-wake formulation in the 
Barnwell and Wahls wall-function approach was only developed for the inner layer to 
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extend the inner layer into the inner portion of the outer layer. Thus, two functions 
characterize the inner layer, whereas only the defect-stream function characterizes the 
outer layer. 

a) Inner-Layer Function: Law-of-the-Wall 

The law-of-the-wall describes the logarithmic behavior of the turbulent velocity 
profile in the inner portion of the boundary layer. This well-known formula is used 
as part of the analytic inner-layer solution in the wall shear stress calculation. Examples 
of the compressible law-of-the-wall formulation and a typical turbulent CFD velocity 
profile are presented in Figure 7. 



Figure 7: Law-of-the-Wall Velocity Profile. 


The region from y + of 30 to 6000 is the logarithmic region and is modeled by the 
law-of-the-wall, equation (76). 
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The form of the law-of-the-wall utilized by Barnwell and Wahls [3] for compressible. 


non-adiabatic flows is: 
Law-of-the-wall: 



The variables « and b are set at 0.41 and 4.9, respectively as suggested by Barnwell [4]. 

b) Inner-Layer Function: Law-of-the-Wake 

The law-of-the-wake has been utilized in this wall-function method for two reasons. 
First of all, the law-of-the-wake allows the effects of the streamwise pressure gradient 
to influence the inner region of the boundary layer. Secondly, the law-of-the-wake 
analytically extends the law-of-the-wall region, thereby increasing the distance between 
the wall and the match point [3]. The law-of-the-wake developed by Barnwell and Wahls 
is only applied to the inner region of the boundary layer and has the following form, 

Law-of-the-wake: 



c) Outer-Layer Function: Defect-Stream Function 

The defect-stream function formulation determines the outer analytic-velocity profile. 
The general definition for the defect-stream function, developed by Clauser [7], in the 
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transformed normal coordinate system, ?/, is: 


d£ = u( 7?) - Up 

di] 


i r 


( 78 ) 


The t) coordinate is defined by the equation, 

1 r Piy) 

^ Jo Poo 


dy 


(79) 


Additionally, the solution developed by Barnwell [3] for the defect-stream function 
at the match point is: 

a l'(? + a ) Af f 1 1 AT \ 

■ r(i + «) 1 (“ rr Nm ) 

(Vm ~ Vw) Ai f 3 ^ ) 

~ —zyr M { a ' 2 ' N ^)} 

The wall-function parameter, a, is defined as, 


df 

di] 


= —Ce~ Nm { _ 

m ( \/{ ( 


(80) 


1 


1 


0 . 1 + 

n 1+2(3 


(81) 


and the transformed coordinate // is: 


V = 


Pw u* 




(82) 


Poo PwU OG&V 

The functions T and M are the Gamma and the Kummer functions, respectively. The 
Gamma function is a generalized factorial function [13]. The Kummer function has the 
following Taylor series expansion, 


M (a, b, N) = 1 + -N + — - 1 ) + 1 H fl + g) 

b b{b+ 1) 2! 6(6+l)(6 + 2) 3! 


(83) 


The other parameters in equation (80) are: 


C = \ Ml a, N w ) + 


1(7 + «) - , f f 1 3 \ 

: n 1 + a) VwM (2 + a ’ 2 ’ Au y 


-1 


® j ) — ' w f ^ 




(84) 
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1 + 2/3 j/u, 


(85) 


e 


JV U , = 


U,’ 


7/u; = 


>w/k 2 

e 

1+2/3 



AU = 


1+2 /3 ( l'}in ij w ) 


u; 


u ./f 


( 86 ) 

(87) 

( 88 ) 


(note: fj w is not r/ evaluated at the wall, but t/ m denotes fj evaluated at the match point). 
These variables are utilized in the defect-stream function to calculate the outer analytic 
velocity at the match point and the match point velocity is then utilized to calculate the 
wall shear stress. 


2.3 Wall Shear Stress 


Once the match point location is determined and the velocity functions are evaluated 
at the match point, the wall shear stress is calculated. The analytic velocities for the 
inner layer and the outer layer are calculated at the match point using the following 
relationships: 

Inner layer formulation: 

“ n, = * [() ( vtn ) + M Dm ) ] (89) 


Outer layer formulation: 


Mm — ^ oo H" M 

ri// 


(90) 


where u* is the shear-stress velocity, g is the law-of-the-wall, h is the law-of-the-wake, 
and ^ is the defect-stream function. The two match point velocity formulations are 
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equated, 


“*b(»m) + A(»+)] = Uoo + U*^- 
and the shear-stress velocity is solved for, yielding, 


u* = 


{•to) + *W) - (fc)v^ U\J 


where the density ratios in the denominator result from the coordinate transformation 
from i] to rj. The shear-stress velocity determines the wall shear stress via the definition, 


T w — p u , U 


thus the wall shear stress is analytically calculated without having to numerically resolve 
the laminar-viscous sub-layer. 


2.4 Slip-Wall Boundary Conditions 

The wall shear stress is calculated analytically using the wall-function method and 
the velocity gradient at the wall is calculated from, 


T w = (m + fit) 


where /q is calculated from Sutherland’s law, equation (21), and fit is calculated from 
the outer-layer turbulence model, the Baldwin-Lomax model, equation (45). Rearranging 


yields, 


du _ Tw 

dy w /'/ + m 
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(95) 



Utilizing the Baldwin-Lomax turbulent eddy-viscosity model and the definition of shear- 


stress velocity, results in the relationship, 

du 

dy 


pw ^ 


(96) 


Pi -|- k Pa C cp F wake 

where the Fueb parameter is approximately 1.0 near the wall. The parameter p s is the 
slip-wall density and is calculated from equation (99). A first order finite-difference 
extrapolation of the velocity gradient is used to calculate the slip-wall velocity. The 
velocity at the first grid point off the wall and the velocity gradient are used to determine 


the slip-wall velocity boundary condition using, 

u s = u(2) - [y ( 2 ) - i/(l)] ^ 


(97) 


where u(2) and y(2) are values at the first grid point off the wall and y(l) is the 
normal coordinate of the wall. The normal-wall velocity is zero, consistent with no fluid 
penetrating the solid surface. Hence, the wall-function velocity boundary conditions for 


the wall are: 


tt(l) = u s ± 0 


(98) 

t’U) = o 

Note that the normal velocity boundary condition on the upper surface (denoted as oo) 
is set to zero (i.e. a streamline) for the test cases examined in this research. 

Along with calculating a slip-wall velocity, a slip-wall density must be calculated 
to be consistent with the slip velocity. The slip-wall density is derivable from Crocco’s 


theorem, yielding, 


P oo 
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(99) 


36 



where the density boundary condition is: 


Pi 1 ) = Ps 


( 100 ) 


The pressure at the wall is determined by satisfying the standard pressure boundary 
condition. 


dp 

0y 


= 0 


( 101 ) 


The numerical boundary condition employed by the finite-difference code is: 


2,, (2) - i,,(3) 

P(l) = 3T= 


( 102 ) 


and was not modified for this research. 

The slip-wall temperature is determined from the slip-wall density, the pressure at 
the wall, and the equation of state as follows, 


T, 


Pi 1 ) 

p a R 


where R is the ideal gas constant, 


R = 287 


J 

k<j ■ K 


1716 


ft- lb 
shig ■ ° R 


Hence, the temperature boundary condition is: 


(103) 


(104) 


T( 1) = T s . (105) 

To account for the streamwise variation in the slip-wall density and its influence on 
the turbulent viscosity, a slip-wall turbulent viscosity is calculated from the wall shear 
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stress, the velocity gradient at the wall, and the laminar viscosity at the wall, in the 


following manner. 




T w 


On 


TTy 

W 


(106) 


where the wall shear stress and the velocity gradient have both been updated. This is 


required to guarantee consistency within the wall-function method. 


2.5 Integration Step Size 

The integration step size is directly proportional to the minimum spacing between 
the grid points in the normal direction. Having increased the minimum spacing between 
two adjacent grid points by implementing the wall-function method, the step size is 
proportionally increased. Some examples of step sizes for both the fully-gridded CFD 
code and the wall-function method are presented in Table 2. 


Table 2: Examples of Non-Dimensional Step Size Values for Turbulent Flat Plate Flow. 



Fully-Gridded CFD 

Wall-Function 

Moo = 2.0 

9 x 10' 8 

3 x 10' 5 

o 

* n 
II 

8 

£ 

1 x 10' 6 

4 x 10 5 

o 

CO 

II 

8 

s 

3 x 10‘ 6 

1 x 10 4 


The Barnwell and Wahls wall-function method has been discussed in detail. The 
match point equation was derived using the Baldwin-Lomax turbulence model consistent 
with the PNS CFD code. The analytic-velocity functions were presented along with the 
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procedure for calculating the wall shear stress and the relevant boundary conditions. The 
effect of the wall-function method implementation on the integrated step size was also 
presented. 
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3. IMPLEMENTATION OF THE 
DEFECT WALL-FUNCTION METHOD 


The prior section discussed in general the solution process for the Barnwell and 
Wahls wall-function method. The current research effort modified the Korte PNS CFD 
code to utilize the Barnwell and Wahls wall-function method. Several modifications to the 
PNS CFD code and the wall-function methodology were implemented to integrate the two 
entities. The procedure for utilizing the wall-function method to generate numerical finite- 
difference solutions is discussed in the next section along with the required modifications 
to the PNS CFD code. 


3.1 Procedure 

The procedure for implementing the wall-function method required that the finite- 
difference CFD code utilize no-slip boundary conditions for the first few solution planes 
to generate initial conditions for the wall-function method, since an initial solution for 
the flowfield must be known to initiate the Barnwell and Wahls wall-function method. 
The initial solution must contain the streamwise velocity profile, u(y), the density profile, 
/>( y), and also specify the parameters, 


Dm 



and u s ■ 


( 107 ) 


For the first computational plane employing the wall-function method, the values for y m 
and u s are unknown, thus the assumed values for these two quantities are the y and u 
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values at the first grid point off the wall, 


f u(2) i = 1 

\ i > 1 


( 108 ) 


y»i , 


r 2/(2) 2 = i 

\ V m 2 -> 1 


( 109 ) 


(where “i” denotes the computational plane being calculated using the wall-function 
method). Subsequent computational planes utilize the slip conditions from the previous 
computational plane. The modified PNS CFD code uses the slip boundary conditions 
in the same manner as the no-slip boundary conditions. Since the wall shear stress has 
been calculated analytically, the first grid point off the wall is moved outward to allow 
for larger grid spacing. 

Several flowcharts were developed to illustrate the steps taken to implement the wall- 
function methodology in the PNS CFD code. The “main program” structure is shown 
in Figure 8 and is discussed in the following text (note: not all the steps/calculations 
are presented in the flowcharts, just the most important ones related to the wall-function 
implementation). 
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Figure 8: The “main program” Flowchart. 


The “main program” flowchart consists of several “call” statements and an “if’ statement 
to either continue marching or to stop. The “call” statements for the predictor integration 
step (labeled “pred”) and corrector integration step (labeled “corr”) are presented in the 
“main program” flowchart as previously discussed in Section 1.4. 

The first subroutine “called” by the “main program” is the “initial” subroutine, shown 
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in Figure 9. 



Figure 9: The “initial” Subroutine Flowchart. 


The “initial” subroutine “reads in” the initial values for several variables, such as the 
stretching factor and the freestream conditions. The spatially uniform initial values for 
all the non-dimensional flowfield variables are set in this subroutine (all cases examined 
in this research were initialized to uniform fields). 

The next subroutine "called" by the “main program” is “code”. This subroutine 
calculates all the flowfield variables from solutions of the PNS equations. The “boundary” 
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subroutine is "called" from the “code” subroutine. The “boundary” subroutine, shown 
in Figure 10, sets the boundary conditions for the flowfield variables, such as velocity, 
density, pressure, and temperature. 
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Note the variable istep and the variable itran are used to initiate the wall-function 
procedure. The variable istep is an the integer number given to each flowfield plane 
and the variable itran is the integer number for the first flowfield plane employing the 
wall-function slip-wall boundary conditions (typically this user specified value has been 
set to 3 in this research). The steps shown on the left side of Figure 10 are for the 
wall-function method implementation and the steps on the right are for the PNS CFD 
code using no-slip boundary conditions. The variables u*(l) and M(l) are set in the 
“wallmain” subroudne for the wall-function method. 
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The next subroutine "called" by the “main program” is “step”. Figure 11. 



Figure 11: The “step” Subroutine Flowchart. 


The variables y* m i n and y* max are the minimum and maximum non-dimensional y values, 
respectively. The first step determines the wall shear stress from either the wall-function 
calculated value or the numerically calculated value. This subroutine calculates the non- 
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dimensional axial step size, dx*, between consecutive flowfield planes and also calculates 
the non-dimensional y-coordinates of the grid at each new plane. 

In this subroutine, the wall-function method calculates the non-dimensional y- 
coordinates slightly different than the original code. The wall-function method deter- 
mines the location of the first grid point off the wall, y*(2), based on a percentage, pm, 
of the match point location, y* m . Typical percentages have ranged from 5 % to 100% for 
this research. The other grid points are calculated using the stretching factor algorithm. 
The reasoning and discussion for the changes in the “step” subroutine are presented in 
Section 3.2. When using the wall-function method, the stretching factor is a constant, 
chosen to keep several grid points in the boundary layer. 

In the “boundary” subroutine, the wall-function method is "called" through the 
subroutine “wallmain”, Figure 12. 



Figure 12: The “wallmain” Subroutine Flowchart. 
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This subroutine calculates the freestream conditions for the speed of sound, velocity, 
laminar viscosity, and density to dimensionalize the flowfield variables before utilizing 
the wall-function method, which was developed for primitive (dimensional) variables. 
The “wallsub” subroutine, shown in Figure 13, is "called” from “wallmain” and is 


basically the Barnwell and Wahls wall-function methodology. 



Figure 13: The “wallsub” Subroutine Flowchart. 
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This subroutine must assume some initial values the very first time wall-function cal- 
culations are performed (i.e. jfirst=0), as discussed earlier in Section 3.1. Subsequent 
calculations use the values from the previous plane (i.e. i-1) as the initial guesses. Sev- 
eral wall-function variables must be calculated using the Barnwell and Wahls equations. 
The variable fj is solved for using the match point equation. After more calculations, 
the velocity gradient at the wall, the slip-wall velocity, and the wall shear stress are de- 
termined. The flowfield variables are returned to “wallmain” to be non-dimensionalized 
and then returned to the “boundary” subroutine. 

The subroutine “fwanal”, shown in Figure 14, is "called" from the “eddy” subroutine, 
which is "called" from the “code” subroutine, to calculate the variables F wa k e and y ma x 
needed by the Baldwin-Lomax turbulence model (note: when referring to the turbulence 
model, y ma x refers to the y value at the maximum F value). Another modification of 
the PNS CFD code required for the implementation of the wall-function methodology 
is to apply the outer-layer turbulent viscosity at all grid points including those in the 
inner-layer [2], This is performed in the “eddy” subroutine. 
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Figure 14: The “fwanal” Subroutine Flowchart. 


This subroutine calculates an analytical grid (labeled “y an ”)- Using this analytical grid. 
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“fwanal” calculates the streamwise analytic-velocity profile using the analytic functions, 
which in turn is used to calculate the vorticity distribution in the boundary layer. The 
vorticity is used to calculate the F(y) function, which determines the F ma x and y ma x 
needed to calculate the F wake parameter. This procedure compensates for the lack of 
resolution in the numerical CFD grid and was found to yield reasonable results. Further 
discussion of the implementation of the “fwanal” subroutine and the relevant changes to 
the wall-function equations is presented in Section 3.3. 


3.2 Defect Wall-Function Gridding Scheme 

After the wall-function method was implemented, the original gridding scheme was 
altered. The numerical grid is developed based on the location of the match point. The 
match point location is calculated in the “wallsub” subroutine and then the first grid point 
off the wall is chosen based on a user specified percentage of the match point location, 
pm, in the “step” subroutine. Moving the first grid point closer to the wall allows for 
more grid points to lie within the boundary layer. Once the first grid point off the wall 
is determined, the remaining grid points are controlled by the stretching factor algorithm. 
Although no comprehensive sensitivity study was performed, typical values for the pm 
parameter and the sfact parameter, that yielded reasonable results, are 0.50 and 1.005, 
respectively. It is important to note that the sfact parameter is critical, since the minimum 
grid spacing controls the integration step size and thus net computational time. 
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3.3 Modified Turbulence Modeling for Defect Wall-Function Method 


Since the numerical grid resolution is reduced in the boundary layer by moving grid 
points farther from the wall, the vorticity for the turbulence model may not be adequately 
specified. To address this issue, an analytical grid is established on which an analytical 
vorticity distribution is calculated. From this data, an analytical F^e and y ma x are 
calculated for the Baldwin-Lomax turbulence model. 

A semi-uniform analytic-grid is generated with the first analytic-grid point off the wall 
located at a y + of 1.0 and a uniform y + analytic-grid spacing from that point outward 
(typical grid spacings have ranged from 5-20 in this research). The velocity for the 
analytic-grid points below the match point is determined from the law-of-the-wall and 
law-of-the-wake. The changes to the relevant wall-function equations are: 
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and the subscript “an” denotes an analytic-grid calculation. The analytical velocity for the 
analytic-grid points above the match point is determined from the defect-stream function 
formulation. The equations are: 
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and pan is the analytically calculated density and is determined from the following 
equation, developed by Barnwell and Wahls [4], as a function of the defect-stream 
function, 

Pan{ytn) = - ' ^ 7 \ - - i ’ ( 117 ) 
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Once the transformation is used, equation (117) becomes a quadratic in p an ■ This 
quadratic is solved and the largest root corresponds to the density profile, (empirically 
verified by comparison with the density profiles generated employing the fully-gridded 
CFD code). Once the analytic-velocity profile is established, the analytic vorticity is 
calculated and the F wakc and y max quantities for the Baldwin-Lomax turbulence model 
are determined. This analytical calculation ot the Baldwin-Lomax parameters allows for 
a more accurate calculation of the outer- layer turbulent eddy viscosity, since there are 
fewer numerical grid points in the boundary layer when implementing the wall-function 
method. 
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The procedure for implementing the Barnwell and Wahls wall-function method has 
been discussed (graphically illustrated in the flowcharts) as well as the modifications to the 
PNS CFD code. The numerical grid generating scheme was altered to take advantage of 
the wall-function implementation allowing larger grid spacing. The analytical calculation 
of the parameters for the outer-layer turbulent eddy-viscosity model was discussed. 
As shall be seen, all of these modifications allowed the PNS CFD code to generate 
engineering accurate results very quickly as compared to the fully-gridded CFD code. 
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4. RESULTS 


The PNS CFD code has been modified to incorporate the Barnwell and Wahls wall- 
function method. The PNS CFD code was validated using laminar flow conditions on a 
flat plate (zero streamwise pressure gradient) to ensure the code was functioning properly 
with respect to theoretical velocity and temperature profiles and by checking conservation 
of mass, momentum, and energy. Several turbulent flat plate (zero streamwise pressure 
gradient) test cases were also investigated, utilizing both the wall-function methodology 
and the fully-gridded methodology (resolved laminar-viscous sub-layer). The resulting 
solutions were compared between these two methods and also to theoretical distributions 
of the wall shear stress and experimental data for the velocity profile to validate both 
methods. Non-zero streamwise pressure gradient cases were also investigated. 


4.1 Laminar Flat Plate Flow (Zero Pressure Gradient) 

a) Velocity and Temperature Profiles 


The PNS CFD code was evaluated with a laminar-viscous formulation to ensure that 
there were no errors in its implementation (without the complicated issues of turbulence 
modeling). The evaluation involved a comparison with laminar boundary-layer profiles 
for flat plate flow developed by Crocco [11], Crocco’s exact solutions for laminar, 
adiabatic flat plate flow use a linear viscosity law (u; = 1) and a Prandtl number of 1.0. 
The parameter u is utilized in the following manner, 
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Crocco’s laminar solutions are plotted in Figures 15 and 16 for several Mach numbers. 


These graphs are obtained from Schlichting, reference [11]. 



Figure 15: Theoretical Velocity Profiles for Laminar Flat Plate Flow. 


T 
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Figure 16: Theoretical Temperature Profiles for Laminar Flat Plate Flow. 


The CFD code was tested with flat plate flow at two different Mach numbers. The 
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domain for the computations is: 

0 < j < 1 
XJ 

0 < y < 1.2 . 

±J 

The boundary conditions for the CFD code are listed below, 

u ( 1 ) = 0 

u{llj ) — f/oo 

u(l) = 0 


( 120 ) 


v(nj) = 0 

T(nj) = T 00 (121) 

P(nj) = p(nj) R T(nj) 
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nj = 121 , 

where the parameter nj is the number of y-grid points in the computational field. The 
pressure boundary condition at the wall is calculated so that the following relation is 
satisfied, 
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The temperature boundary conditions for the wall are: 


Adiabatic wall: 


( 122 ) 


Pseudo-adiabatic wall: 
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Pseudo-adiabatic denotes that the wall temperature is set to be the adiabatic wall temper- 
ature computed by equation (72), whereas adiabatic denotes that the temperature gradient 
at the wall is explicitly zero. Cases A- 1 and A- 2 have a freestream Mach number of 2.0, 
while cases B-l and B-2 have a freestream Mach number of 5.0. Cases A-l and B-l 
used an adiabatic wall temperature boundary condition and cases A-2 and B-2 used a 
pseudo-adiabatic wall temperature boundary condition. The other inflow conditions are 


listed in Table 3. 

Table 3: Laminar Flat Plate Inflow Conditions. 


Parameter 

Case A-l 

Case A-2 

Case B-l 

Case B-2 

Mqc 

2.0 

2.0 

5.0 

5.0 

Re L ! 

1.5 x 10 6 

1.5 x 10 6 

1.5 x 10 6 

1.5 x 10 6 

Too 

222.0 K 

222.0 K 

222.0 K 

222.0 K 

T w 

399.6 K 

399.6 K 

1332.0 K 

1332.0 K 

Pr 

1.0 

1.0 

1.0 

1.0 

UJ 

1.0 

1.0 

1.0 

1.0 

L (length) 

1.0 m 

1.0 m 

1.0 m 

1.0 m 

Temp B.C. 

adiabatic 

pseudo-adiab. 

adiabatic 

pseudo-adiab. 

Temp B.C. 

dT _ n 

T{\) - T aw 

O 

II 

T(l) = T aw 
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The velocity profiles for cases A-l and A-2 are presented in Figure 17. 



Figure 17: Laminar, Velocity Profiles (Case A). 


Three velocity curves are presented in Figure 17. One curve represents Crocco theory. 
The other two curves are calculated using the PNS CFD code. One uses an adiabatic 
wall temperature boundary condition, whereas the other uses a pseudo-adiabatic wall 
temperature boundary condition. The CFD code uses equation (125) to numerically 
calculate the non-dimensional adiabatic wall temperature boundary condition, T*(l), 


n i) 


T*( 2) + 


[T*(2) -T*(3)] 


»♦( 3)-y»(2)' 


Jf*(2)-r(l ). 
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y(3)-y(2> 


)} 


(125) 


This equation numerically approximates equation (123). The third curve uses a pseudo- 
adiabatic wall temperature boundary condition, equation (124). The two CFD curves 
compare well with the theory of Crocco. The maximum percent error between the CFD 
curves and the theory occurs around ?/= 5 and is only 2.3%. Percent error is defined in 
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the following manner. 


%Err = M//,eorg - l ^££2 

u theory 


( 126 ) 


The temperature profiles for cases A-l and A-2 are presented in Figure 18. 



Figure 18: Laminar, Temperature Profiles (Case A). 


The adiabatic boundary condition does not predict the correct adiabatic wall temperature, 
calculated using equation (72). For case A, the non-dimensional adiabatic wall temper- 

'T' 

ature, is 1.8. The pseudo-adiabatic boundary condition closely approximates the 
adiabatic wall temperature boundary condition given by equation (123) by producing a 
near zero temperature gradient at the wall. The use of the numerical adiabatic boundary 
condition, equation (125), results in a percent difference of 0.03% between the wall tem- 
perature (point 1) and the first point off the wall (point 2). Using the pseudo-adiabatic 
boundary condition also gives a percent difference of 0.03% between the two correspond- 
ing points. This implies that the pseudo-adiabatic boundary condition approximates the 
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adiabatic boundary condition. Percent difference is defined as. 


%Diff 


T{l)-T{2) 

W) 


( 127 ) 


The two CFD curves closely approximate Crocco theory. The maximum percent error 
between the CFD curves and the theory occurs around rj=6 and is 3.7%. 

The velocity profiles for cases B-l and B-2 are presented in Figure 19. 



Figure 19: Laminar, Velocity Profiles (Case B). 

The CFD code closely models Crocco theory. The percent error between the pseudo- 
adiabatic CFD curve and the theory is a maximum around i]= 15 and is 3.2%. There is 
more separation between the adiabatic and pseudo-adiabatic curves between rj= 13 and 
r/=16 in case B than for case A, but both curves are close to Crocco theory. 
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The temperature profiles for cases B-l and B-2 are presented in Figure 20. 



Figure 20: Laminar, Temperature Profiles (Case B). 

Using the numerical adiabatic boundary condition, equation (125), results in a percent 
difference of 0.2% between the wall temperature (point 1) and the first point off the wall 
(point 2). Using the pseudo-adiabatic boundary condition also gives a percent difference 
of 0.2% between corresponding points, implying an adequate adiabatic boundary con- 
dition. The pseudo-adiabatic boundary condition models Crocco theory closer than the 
adiabatic boundary condition from ?/=0 to r/=13. The maximum percent error between 
the pseudo- adiabatic curve and the theory occurs around 77=15 and is 32.0%. This latter 
error could be reduced via grid resolution, but this effort was deemed unimportant for 
this study. 
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b) Conservation Laws 


To ensure that the code is functioning properly with respect to conservation of mass, 
momentum, and energy, these quantities were investigated for the laminar cases. The 
conservation of mass is given in the following equation, the mass flow rate per unit width, 

y max 

p{y)u(y) dy , (128) 

where m is the mass flow rate, w is the width, and y ma x is the y-location of the outer 
boundary (which is a stream surface). For the two-dimensional cases dealt with in this 
research, w is 1.0. The non-dimensional mass flow rate is given in the following equation, 



where the * represents a non-dimensional quantity and L is the characteristic length of 
the flat plate. The non-dimensional mass flow rate was calculated and found to be a 
constant in all four cases and approximately equal to 1.2 across the entire length of the 
flat plate, implying that mass is conserved. The non-dimensional mass flow rate for the 
entire flowfield and the boundary layer at the trailing edge (i.e. x=1.0 m) are shown 
in the following table, as well as the percentage of the entire flowfield non-dimensional 



mass flow rate in the boundary layer. 

Table 4: Non-Dimensional Mass Flow Rates at Trailing Edge. 


Test Case 

Entire Flowfield 

Boundary Layer 

% of Entire Field 

Case A- 1 

1.200 

0.085 

7.08% 

Case A-2 

1.200 

0.085 

7.08% 

Case B-l 

1.200 

0.233 

19.4% 

Case B-2 

1.200 

0.235 

19.6% 
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The percent error between the non-dimensional mass flow rate at the leading edge (x=0.0 
m) and the trailing edge is presented in Table 5 for all cases. Percent error is defined as. 


%Err = 


m ; =1 - m;_ 


x=0 


m 


x=\ 


(130) 


Table 5: Percent Errors for Non-Dimensional Mass Flow Rate. 


Test Case 

Percent Error for O.Oc-^cl.O 

Case A-l 

0.001% 

Case A- 2 

0.001% 

Case B-l 

0.001% 

Case B-2 

0.001% 


The small percent errors imply that the non-dimensional mass flow rate is constant, 
thus mass is conserved. 

The conservation of momentum using the stream thrust approach is presented below. 
Momentum is conserved if the following condition is met, 


AF * = Drag 


(131) 
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The stream thrust per unit width is: 
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where F is the stream thrust. The non-dimensional stream thrust is given below. 


F* 
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The non-dimensional stream thrust and drag for cases A and B are given in Figures 21, 
22, 23, and 24. 



Figure 21: Conservation of Momentum (Case A-l). 
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Figure 24: Conservation of Momentum (Case B-2). 


The change in stream thrust between two streamwise locations and the drag on the surface 
between those streamwise locations are approximately equal for all four cases and differ 
by less than 7% across the entire length of the flat plate, implying that momentum is 
conserved. The maximum percent error between the curves over the length of the flat 
plate and the percent error at the trailing edge for all cases are given in Table 6. 


Table 6: Percent Errors for Non-Dimensional Stream Thrust Approach. 


Test Case 

Maximum % Error 


Case A-l 

6.49% at f =0.06 

5.52% 

Case A- 2 

6.45% at f =0.06 

5.52% 

Case B-l 

6.52% at £=0.12 

4.45% 

Case B-2 


4.21% 


The conservation of momentum is also checked using the momentum integral equa- 
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tion for two-dimensional compressible flow over a flat plate as. 


Afl _ Cf _ Ty 
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where 6 is the momentum thickness. 
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and Cf is the coefficient of friction. The non-dimensional momentum integral equation is: 


A 0_ 
Ax 


* 



(138) 


The results from the non-dimensional momentum integral equation for cases A and B 
are given in Figures 25, 26, 27, and 28. 



Figure 25: Conservation of Momentum (Case A-l). 
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Figure 26: Conservation of Momentum (Case A-2). 



Figure 27: Conservation of Momentum (Case B-l). 
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Figure 28: Conservation of Momentum (Case B-2). 


The change in momentum thickness between two streamwise locations and the non- 
dimensional wall shear stress on the surface between those streamwise locations is 
presented in the four graphs and differs by less than 10% across the flat plate for j>0.06 
for case A and for j>0.18 for case B, implying that momentum is conserved in these 
regions. The maximum percent error between the curves over the flat plate and the 
percent error at the trailing edge for all cases are given in Table 7. 


Table 7: Percent Errors for Momentum Integral Equation Approach. 


Test Case 

Maximum % Error 

uzEsnai 

Case A-l 


0.19% 

Case A-2 

0.62% at f =0.21 

0.18% 

Case B-l 


2.15% 

Case B-2 

9.77% at f =0.69 

2.18% 
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The conservation of energy is given in the following equation, the product of the 


mass flow rate and the total enthalpy per unit width. 
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where H t is the total enthalpy. This product is denoted as the energy flux throughout the 
rest of the text. The non-dimensional energy flux is given as. 
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where 


C* 
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(7"l)^i 
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Energy is conserved if the energy flux is a constant (for an adiabatic wall case). The 
non-dimensional energy flux was calculated and found to be a constant approximately 
equal to 1.35 for case A and approximately equal to 0.72 for case B across the entire 
length of the flat plate, implying that energy is conserved. The non-dimensional energy 
flux for the entire flowfield and the boundary layer at the trailing edge are shown in the 
following table, as well as the percentage of the entire flowfield non-dimensional energy 
flux in the boundary layer. 


Table 8: Non-Dimensional Energy Fluxes at Trailing Edge. 


Test Case 

Entire Flowfield 

Boundary Layer 

% of Entire Field 

Case A-l 

1.350 

0.191 

14.1% 

Case A-2 

1.350 

0.193 

14.3% 

Case B-l 

0.720 

0.281 

39.0% 

Case B-2 

0.720 

0.283 

39.3% 
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The percent error between the non-dimensional energy flux at the leading edge and the 
trailing edge is presented in Table 9 for all cases. 


Table 9: Percent Errors for Non-Dimensional Energy Flux. 


Test Case 

Percent Error for 0.0<j<1.0 

Case A- 1 

0.001% 

Case A-2 

0.001% 

Case B-l 

0.0003% 

Case B-2 

0.0001% 


The small percent errors in the non-dimensional energy flux imply that it is constant and 
energy is conserved. The PNS CFD code has been shown to be operating correctly for the 
laminar cases by comparison with theoretical profiles and by checking for conservation 
of mass, momentum, and energy. 


4.2 TUrbulent Flat Plate Flow (Zero Pressure Gradient) 

The wall-function method was implemented into the CFD code and a flat plate 
model was used to test the wall-function method versus the fully-gridded CFD code, 
an Eckert reference method, and specific experimental data. The fully-gridded CFD test 
cases have resolved laminar-viscous sub-layers and no-slip boundary conditions. The 
pseudo-adiabatic wall temperature boundary condition, equation (124), is used for both 
the wall-function cases and the fully-gridded CFD cases. The pseudo-adiabatic test case 
inflow conditions are listed in Table 10, (where W-F denotes wall-function cases and 
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F-G denotes fully-gridded CFD cases). 


Table 10: Turbulent Flat Plate Pseudo-Adiabatic Inflow Conditions. 


Parameter 

Case C 

Case D 

Case E 

Mqo 

2.0 

5.0 

8.0 

ReL 

20.0 x 10 6 

15.0 x 10 6 

20.0 x 10 6 

Too 

222.0 K 

100.0 K 

150.0 K 

T w 

381.2 K 

548.1 K 

1870.9 K 

Pr 

0.72 

0.72 

0.72 

Pr, 

0.9 

0.9 

0.9 

L (length) 

1.0 m 

1.0 m 

1.0 m 

sfact (W-F) 

1.001 

1.001 

1.002 

sfact (F-G) 

varies with x 

varies with x 

varies with x 

CFL (W-F) 

0.25 

0.4 

0.5 

CFL (F-G) 

0.1 

0.1 

0.1 

pm 

0.25 

0.5 

0.7 


The domain of the computations is the same as shown in equation (120). The boundary 
conditions are the same as those in equations (121), (122) and (124) for the fully-gridded 
CFD case, with the changes listed in Section 2.4 for the wall-function method. 


a) Conservation Laws 

Similar to the laminar test case, the conservation laws were checked to ensure that 
the PNS CFD code with no-slip boundary conditions and with slip boundary conditions 
is functioning properly with respect to conservation of mass, momentum, and energy. 
The non-dimensional mass flow rate for all three cases was calculated and found to be 
a constant around 1.2 for the entire length of the flat plate for both the wall-function 
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and fully-gridded cases, implying that mass is conserved. The non-dimensional mass 


flow rate for the entire flowfield and the boundary layer at the trailing edge are shown 
in the following table, as well as the percentage of the non-dimensional mass flow rate 
in the boundary layer. 


Table 11: Non-Dimensional Mass Flow Rates at the Trailing Edge. 


Test Case 

Entire Flowfield 

Boundary Layer 

% in B. L. 

Case C (W-F) 

1.200 

0.0186 

1.55% 

Case C (F-G) 

1.200 

0.0235 

1.96% 

Case D (W-F) 

1.200 

0.0167 

1.39% 

Case D (F-G) 

1.200 

0.0081 

0.68% 

Case E (W-F) 

1.200 

0.0173 

1.44% 

Case E (F-G) 

1.200 

0.0084 

0.70% 


The percent error between the non-dimensional mass flow rate at the leading edge and 
the trailing edge is presented in Table 12. 


Table 12: Percent Errors for Non-Dimensional Mass Flow Rate. 


Test Case 

Percent Error for O.Oc^cl.O 

Case C (W-F) 

0.029% 

Case C (F-G) 

0.015% 

Case D (W-F) 

0.016% 

Case D (F-G) 

0.002% 

Case E (W-F) 

0.007% 

Case E (F-G) 

0.001% 
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The conservation of momentum using the stream thrust approach for cases C, D, and 


E are presented in the following figures. 



Figure 29: Conservation of Momentum (Case C). 



Figure 30: Conservation of Momentum (Case D). 
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Figure 31: Conservation of Momentum (Case E). 


Momentum is conserved if equation (131) is satisfied. For all three wall-function cases, 
the two curves differ by less than 50% for the entire length of the flat plate. The distance 
between streamwise locations (for the integration of the drag and the change in stream 
thrust) differs for the fully-gridded and wall-function cases, implying that the curves for 
the fully-gridded case should be consistent with each other, but should not be consistent 
with the curves of the wall-function method. For the fully-gridded CFD case, the curves 
match well for case D and E, but diverge for case C. The maximum percent error between 
the curves over the length of the flat plate and the percent error at the trailing edge for 
all cases are given in Table 13. 
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Table 13: Percent Errors for Non-Dimensional Stream Thrust Approach. 


Test Case 

Maximum % Error 

BIBHI 

Case C (W-F) 

49% at f =1.0 

49% 

Case C (F-G) 

28% at £=1.0 

28% 

Case D (W-F) 

bbhh 

42% 

Case D (F-G) 

6% at £=1.0 

6% 

Case E (W-F) 

40% at £=0.84 

35% 

Case E (F-G) 

5% at £=1.0 

5% 


Another momentum check for the flat plate case is the momentum integral equation 


given in equation (136). The results are plotted in the following figures. 



Figure 32: Conservation of Momentum (Case C). 
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Figure 33: Conservation of Momentum (Case D). 



Figure 34: Conservation of Momentum (Case E). 


Equation (138) must be satisfied for momentum to be conserved. For the wall-function 
cases, the two curves match closely, although the curve for the term is somewhat 
erratic, especially for case E where there is a spike in the curve near ^=0.77. The 
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fully-gridded CFD curves match well, except for case C. The maximum percent error 
between the curves over the flat plate, for j- >0.2, and the percent error at the trailing 
edge for all cases are given in Table 14. It should be noted that the spatial momentum 
variation is the difference between two large numbers and hence is difficult to accurately 
predict. 


Table 14: Percent Errors for Momentum Integral Equation Approach. 


Test Case 

Maximum % Error 

% Error at x=1.0 

Case C (W-F) 

20% at f = 0.22 

7% 

Case C (F-G) 

271% at f= 0.25 

61% 

Case D (W-F) 

23% at f = 0.94 

13% 

Case D (F-G) 

34 % at x=0.72 

29% 

Case E (W-F) 

85% at ^=0.77 

49% 

Case E (F-G) 


27% 


Energy is conserved if the non-dimensional energy flux, equation (140), is constant. 
The non-dimensional energy flux is nearly constant for both the fully-gridded CFD case 
and wall-function method, implying that energy is conserved. The non-dimensional 
energy flux for the entire flowfield, the non-dimensional boundary layer energy flux at 
the trailing edge, and the percentage of the non-dimensional energy flux in the boundary 
layer are shown in the following table. 
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Table 15: Non-Dimensional Energy Fluxes at the Trailing Edge. 


Test Case 

Entire Flowfield 

Boundary Layer 

% in B. L. 

Case C (W-F) 

1.350 

0.0417 

3.09% 

Case C (F-G) 

1.350 

0.0529 

3.92% 

Case D (W-F) 

0.720 

0.0200 

2.78% 

Case D (F-G) 

0.720 

0.0097 

1.35% 

Case E (W-F) 

0.647 

0.0187 

2.89% 

Case E (F-G) 

0.647 

0.0091 

1.41% 


The percent error between the energy flux at the leading edge and the trailing edge is 
given in Table 16. 


Table 16: Percent Errors for Non-Dimensional Energy Flux. 


Test Case 

Percent Error for 0.0<j<1.0 

Case C (W-F) 

0.032% 

Case C (F-G) 

0.015% 

Case D (W-F) 

0.019% 

Case D (F-G) 

0.002% 

Case E (W-F) 

0.009% 

Case E (F-G) 

0.001% 


In summary, these turbulent calculations provide a verification of the PNS CFD code’s 
ability to preserve, with reasonable accuracy, the flux related quantities for both the 
fully-gridded CFD case and the wall-function method. 


b) Wall Shear Stress 

The analytically calculated wall shear stress from the wall-function methodology 
was compared to the wall shear stress calculated with an Eckert reference method and 
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the wall shear stress calculated with the fully-gridded CFD code. An Eckert reference 
method [14] was utilized to calculate wall shear stress distribution for flat plate flow. 
The inputs needed for the code are listed below. 

Moo , P 0 , T 0 , T w , and x , (142) 


where the subscript “o” denotes a stagnation condition and x is the streamwise distance 
from the leading edge of the flat plate. These quantities are used to calculate the 
freestream conditions of temperature, pressure, density, speed of sound, and velocity 
using the following set of equations, 


T = T 

■*- OC — 1 O 


Poo — 1*0 


1 + 


i; 1 

9 


+ M 


-(^r) 


poo — 


Poo 

1FK 


= n /7 


f'oc — -1'OO^OC 1 


(143) 

(144) 

(145) 

(146) 

(147) 


where R is the ideal gas constant and aoo is th e freestream speed of sound. The freestream 
viscosity is determined from Sutherland’s law of the form. 


3 


[Tool 

2 

' T'i + 5 ‘ 

Tv J 


1 oo + S _ 


(148) 


where 

tr k (j 

/<i = 1.7891 x 10 — — 

7)i sec 

1\ = 288.16 I\ 


(149) 


S = 110.4 K 
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The recovery temperature for turbulent flow is: 


T t = + 


( 150 ) 


and the reference temperature is calculated using Eckert’s formula. 


T re/ = 0.28 + 0.50 T w + 0.22 T r 


(151) 


The reference density is determined using the equation of state. 


H T ref 


(152) 


and the reference viscosity is calculated from Sutherland’s law, equation (148), 


f l rcf — /M 


Ire /] 2 \ T\ + S 


Tl J [Tref + 


(153) 


Employing the reference quantities, the skin friction coefficient for turbulent flow is 
calculated from the following equation, 


0.001 \prefY T/*re/ 


(Re,)* IP 


(154) 


where 


PociRoo-e 


(155) 


Wall shear stress is calculated from the skin friction coefficient utilizing the relationship, 


T W — .y C f Poo R OO 


(156) 
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The comparison between the wall shear stresses of the wall-function method, the fully- 
gridded CFD case, and the Eckert reference method are presented in the Figures 35-37. 



Figure 35: Wall Shear Stress Distribution (Case C). 



Figure 36: Wall Shear Stress Distribution (Case D). 
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Figure 37: Wall Shear Stress Distribution (Case E). 


The wall-function method matched reasonably well with both the Eckert reference method 
and the fully-gridded CFD case, especially for case C, with a percent difference between 
the wall-function case and the fully-gridded CFD case at the trailing edge of only 15%. 
Cases D and E did not match as well, but are adequate approximations with percent 
errors of 25% and 40%, respectively at the trailing edge. 


c) Fully-Gridded CFD Comparison 

The velocity profiles generated with the wall-function methodology are compared (at 
the same x-location) to those generated by the fully-gridded CFD code in the following 
figures. The velocity profiles for the entire field and the velocity profiles for just the 
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boundary layer are presented 



Figure 38: Velocity Profiles for Entire Field (Case C). 



Figure 39: Velocity Profiles for Entire Field (Case D). 
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Figure 40: Velocity Profiles for Entire Field (Case E). 


The velocity profiles for the entire field are shown to illustrate the large gradients defining 
the boundary layer and to emphasize the spatial scale of the boundary layer (i.e. a majority 
of the flowfield is in the freestream). For clarity, the boundary-layer velocity profiles are 
presented in the following figures and demonstrate the Barnwell and Wahls wall-function 
method’s ability to capture the structure of the boundaiy-layer velocity field. 
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Figure 42: Velocity Profiles (Case D). 
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Figure 43: Velocity Profiles (Case E). 


The velocity profiles for the boundary layer for all cases match well. The wall-function 
velocity profile for case C differs from the fully-gridded CFD case in the outer part of 
the boundary layer, but by less than 3%. The profiles generated by the wall-function 
method for cases D and E match with the fully-gridded CFD case everywhere except for 
the wall point (which is not supposed to match). 

The temperature profiles in the boundary layer are also compared between the wall- 
function method and the fully-gridded CFD case, again at the same x-location. 
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Figure 44: Temperature Profiles (Case C). 



Figure 45: Temperature Profiles (Case D). 
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Figure 46: Temperature Profiles (Case E). 


The results for the temperature profiles in the boundary layer are similar to the velocity 


profiles. All three cases match very well throughout the entire boundary layer. 


The boundary-layer density profiles are also compared at the same x-location in the 
following figures. 
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Figure 49: Density Profiles (Case E). 


Again, all three wall-function cases match well in the inner part of the boundary layer 
with the fully-gridded CFD cases, but not as well in the outer part of the boundary layer. 
The implemented wall-function methodology generates profiles of velocity, temperature, 
and density that compare with those obtained from the fully-gridded PNS CFD code. 

To ensure that the modified PNS CFD code is able to yield consistent solutions (as 
compared to those generated by the fully-gridded CFD code), results from cases C, D, and 
E are presented detailing the pressure contours and indirectly, the shock wave angles. 
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Figure 50: Pressure Contours (Case C, CFD). 
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Figure 51: Pressure Contours (Case C, Wall-Function). 
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Figure 52: Pressure Contours (Case D, CFD). 



Figure 53: Pressure Contours (Case D, Wall-Function). 
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Figure 55: Pressure Contours (Case E, Wall-Function). 


These results indicate that the wall-function method reproduces similar pressure contours 
and shock-wave angles (as compared to the fully-gridded CFD case). 
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d) Experimental Data Comparison and Computational Time Required 

The experimental data collected by Mabey et al. [15] was used as a test case to 
analyze the velocity profile produced by the wall-function method and the fully-gridded 
CFD case. The freestream conditions for the experimental case are listed in Table 17. 


Table 17: Experimental Freestream Conditions. 


Parameter 

Experimental Conditions 

Moo 

4.5 

ReL 

27.9 x 10 6 

Too 

62.8 K 

T w 

292.4 K 

Pr 

0.72 

Pr, 

0.9 

L (length) 

1.0 m 

sfact (W-F) 

1.001 

sfact (F-G) 

varies with x 

CFL (W-F) 

0.25 

CFL (F-G) 

0.1 

pm 

0.25 
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The boundary-layer velocity profiles (at the same x-location) are presented in Figure 56. 



Figure 56: Velocity Profiles for the Experimental Test Case. 


Both the wall-function velocity profile and the fully-gridded CFD velocity profile are a 
reasonable approximation of the experimental data, and lend credibility to both techniques. 

The main reason for incorporating the Barnwell and Wahls wall-function methodology 
into the PNS CFD code is to increase computational speed. The computational time 
required to generate solutions from j=0 to j=l on the Sabre computer, a Cray Y- 
MP resident at NASA Langley Research Center, for the fully-gridded and wall-function 
cases are tabulated below. 


Table 18: Computational Time Required. 


Case 

Fully-Gridded CFD 

Wall-Function 

C 

16729 seconds 

554 seconds 

D 

1330 seconds 

141 seconds 

E 

605 seconds 

1 15 seconds 


97 
















The wall-function method utilizing an identical number of grid points as the fully-gridded 
PNS CFD code decreased the computational time required by a factor of 30 for the Mach- 
2 case, by a factor of 9 for the Mach-5 case, and by a factor of 5 for the Mach-8 case. 

e) Implementation Issues 

An analysis to address the significance of a small numerical change in the mass flow 
rate near the wall (in the stream wise direction) due to the non-physical slip-wall velocity 
boundary condition, upon implementation of the wall-function method (i.e. the slip-wall 
velocity and density), utilizes a procedure similar to the displacement thickness method 
for viscous/inviscid interactions. The displacement thickness, 8 * , physically represents 
the distance a wall must be moved into the flow for an inviscid flow analysis to accurately 
represent the retarded mass flow of a viscous problem [ 12). Conservation of mass requires 
that the stream function at the edge of the boundary layer be the same in both the viscous 
and inviscid cases. To illustrate the displacement thickness for an arbitrary boundary 
layer profile, Figure 57 is presented. 
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Figure 57: Illustration of Boundary- Layer Displacement Thickness. 


The condition for which shaded area A is equal to shaded area B defines the displacement 
thickness. The governing equations for the displacement thickness are: 


l 


ye 



Pc Ue dy 


(157) 


and 


6 * 


y* 


p U 


dy 


(158) 


J U V pc«t/ 

where “e” denotes a quantity at the edge of the boundary layer and 8* is the displacement 
thickness characterizing the mass How rate deficiency in the boundary layer [6]. This 
process allows viscous flow to modeled using an inviscid analysis. 

A similar procedure is utilized to address the wall local mass flow rate addition 
resulting from the implementation of the wall-function method and is graphically repre- 
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sented in Figure 58. 



Figure 58: Illustration of Wall-Function Boundary-Layer Displacement Thickness. 


The analytical profile is assumed to be a representation of the exact solution to the flow- 
field from the wall out to the match point and the wall-function method profile is assumed 
to be a representation of the exact solution from the match point outwards. Conserva- 
tion of mass requires that the stream function at the match point be the same as for the 
exact solution. Thus, the wall-function method computation needs to be initiated from 
a displaced wall (i.e. y(l)=yo), where yo is a type of displacement thickness computed 
by requiring the shaded areas A and B in Figure 58 to be equal, or equivalently the 
following equation is solved, 


[Dm fVm 

/ Pan «<wi t hj = / Pm d'J 

J 0 J a o 


(159) 


where “an” denotes the analytical profile and m denotes the match point. The variable 


yo is solved for via an iterative process and the resulting velocity profile (for case C) 
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numerically generated employing yo as the wall location is presented in Figure 59. 



Figure 59: Velocity Profiles (Case C). 


This modification to the Barnwell and Wahls methodology yielded a reduction in the slip- 
wall velocity of approximately 13%. It is interesting to note that the resulting velocity 
profiles are slightly degraded as compared to those generated with the original Barnwell 
and Wahls method However, the corresponding wall shear stress distributions (presented 
in Figure 60) compare favorably with the fully-gridded CFD method 
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Figure 60: Wall Shear Stress Distribution (Case C). 


This effort did not yield significant improvements to the original method and was not 
pursued further due to lack of time. 


4.3 Turbulent Corner Flow (Non-Zero Pressure Gradient) 

a) Expansion Corner 

To demonstrate the applicability of the Barnwell and Wahls method for flowfields 
with streamwise pressure gradients, an expansion corner flowfield is examined. The 
flowfield starts on a flat plate (zero pressure gradient) and is integrated out to the -£=0.2 
location, where a 2.5° downward sloping ramp is located. This test case is evaluated 
with the freestream conditions of case C and is an approximation to the edge conditions. 
The fully-gridded CFD case and wall-function method for the expansion comer case are 
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presented in Figure 61, as well as the flat plate cases. 



o I 1 I J 1 1 1 1 i i 1 1 1 t 1 I 1 I » * ' 

0.00 0.25 0.50 0.75 1.00 

x [m] 


Figure 61: Wall Shear Stress Distribution for the Expansion Corner. 

The wall shear stress distribution for the expansion corner cases decreases rapidly at the 
expansion corner (as expected). The wall-function method expansion corner case has a 
similar trend as the fully-gridded CFD case. The pressure contours for the expansion 
corner wall-function case are presented in Figure 62 and clearly illustrate the leading 
edge compression field as well as the expansion field generated at the corner. 
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Figure 62: Pressure Contours for Expansion Corner. 


b) Compression Corner 

The wall-function method is also applicable to compression corner flowfields. to 
demonstrate this, a flowfield is tested on a flat plate (zero pressure gradient) out to the 
■£=0.2 location and then a 2.5° upward sloping ramp is encountered. The freestream 
conditions are again taken from case C. The fully-gridded CFD case and wall-function 
method for the compression corner case are presented in Figure 63, as well as the flat 
plate cases. 
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Figure 63: Wall Shear Stress Distribution for the Compression Corner. 


As expected, the wall shear stress distributions for the compression corner cases increase 
at the comer and both methods yield similar trends The pressure contours for the 
compression corner wall-function case is presented in Figure 64 and illustrates two 
compression systems: (1) at the leading edge and (2) at the corner. 
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Pressure Contours 
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Figure 64: Pressure Contours for Compression Comer. 


The PNS CFD code has been validated for both the laminar and turbulent flat plate 
cases with respect to conservation and theoretical profiles (Crocco for the former and 
Eckert for the latter) as well as the experimental data of Mabey for the turbulent case. 
The wall-function method produced numerical solutions quickly to within engineering 
accuracy as compared with the fully-gridded CFD code as well as generating solutions 
that compare with theoretical Eckert distributions and the experimental data of Mabey. 
The wall-function method has also been shown to function for non-zero streamwise 
pressure gradient cases. 
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5. SUMMARY 


The original PNS CFD code developed by Korte was modified to utilize the Barnwell 
and Wahls wall-function methodology for turbulent flowfield analysis. The original 
gridding scheme of the PNS CFD code was modified to account for the match point 
location and to consistently adapt the numerical grid point locations based upon the 
location of the first grid point off the wall. Consequently, the reduction in boundary- 
layer grid resolution required a direct analytical calculation of the vorticity (utilized by 
the turbulence model). Additionally, the Barnwell and Wahls theory was modified to 
incorporate the Baldwin-Lomax turbulence model consistent with the PNS CFD code 
(note: this research is the first to modify this wall-function theory for an outer-layer 
turbulence model other than the Clauser model). 

The PNS CFD code was validated with laminar-viscous boundary layer (Mach-2 
and Mach-5) solutions, to demonstrate that the code was functioning properly (without 
turbulence modeling), since the algorithm was shown to generate solutions consistent 
with the laminar theory of Crocco, as well as conserving the mass, momentum, and 
energy fluxes. 

Also, the PNS CFD code was validated for turbulent flow utilizing the Baldwin- 
Lomax turbulence model. A flat plate test case at Mach-4.5 was utilized to test the 
fully-gridded CFD code (resolved laminar-viscous sub-layer). The fully-gridded CFD 
code produced data closely matching the experimental Mach-4.5 data of Mabey (at the 
same x-location). The fully-gridded CFD code also conserved mass, momentum, and 
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energy (to within engineering accuracy) as well as approximately reproducing the wall 
shear stress distributions derived from an Eckert reference method in the Mach number 
range of 2.0 to 8.0. 

Results from the modified code (wall-function methodology) were investigated for 
three turbulent adiabatic flat plate (zero streamwise pressure gradient) test cases and 
two (non-zero streamwise pressure gradient) corner flow test cases. The test cases 
utilized Mach numbers between 2.0 and 8.0, as well as Reynolds numbers between 
15 million and 20 million (per meter). The wall-function method analytically generated 
wall shear stress distributions which were in reasonable agreement with those generated 
from the fully-gridded CFD code and an Eckert reference method. The implemented 
wall-function methodology also generated boundary-layer profiles consistent with the 
fully-gridded CFD code and the experimental data of Mabey. It is important to note that 
an order of magnitude increase in computational speed was obtained employing the wall- 
function method (as compared to the fully-gridded CFD code) and conservation of mass, 
momentum, and energy was shown to be adequate. The two-dimensional spatial pressure 
contours (flat plate test cases) matched closely between the fully-gridded CFD cases and 
those generated from the wall-function method as well as the pressure contours for the 
non-zero streamwise pressure gradient cases. Overall, the resultant trends in wall shear 
stress distributions from the non-zero pressure gradient cases were shown to approximate 
those of the fully-gridded CFD code. 

The resulting modified code, with the implemented wall-function model, is envi- 
sioned to be an engineering tool applicable to complex turbulent aerodynamic design 
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studies based on its enhanced capability to yield engineering data at a vastly reduced 
computational cost. Configurations can be quickly analyzed, and then more stringent 
finite-difference solution methods, using highly resolved grids, can be subsequently ap- 
plied to validate the engineering design. 
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